home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / gnu / glibc108.gz / glibc108 / glibc-1.08.1 / sysdeps / ieee754 / __drem.c < prev    next >
C/C++ Source or Header  |  1993-11-20  |  3KB  |  106 lines
  1. /* Copyright (C) 1992 Free Software Foundation, Inc.
  2. This file is part of the GNU C Library.
  3.  
  4. The GNU C Library is free software; you can redistribute it and/or
  5. modify it under the terms of the GNU Library General Public License as
  6. published by the Free Software Foundation; either version 2 of the
  7. License, or (at your option) any later version.
  8.  
  9. The GNU C Library is distributed in the hope that it will be useful,
  10. but WITHOUT ANY WARRANTY; without even the implied warranty of
  11. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12. Library General Public License for more details.
  13.  
  14. You should have received a copy of the GNU Library General Public
  15. License along with the GNU C Library; see the file COPYING.LIB.  If
  16. not, write to the Free Software Foundation, Inc., 675 Mass Ave,
  17. Cambridge, MA 02139, USA.  */
  18.  
  19. /*
  20.  * Copyright (c) 1985 Regents of the University of California.
  21.  * All rights reserved.
  22.  *
  23.  * Redistribution and use in source and binary forms are permitted provided
  24.  * that: (1) source distributions retain this entire copyright notice and
  25.  * comment, and (2) distributions including binaries display the following
  26.  * acknowledgement:  ``This product includes software developed by the
  27.  * University of California, Berkeley and its contributors'' in the
  28.  * documentation or other materials provided with the distribution and in
  29.  * all advertising materials mentioning features or use of this software.
  30.  * Neither the name of the University nor the names of its contributors may
  31.  * be used to endorse or promote products derived from this software without
  32.  * specific prior written permission.
  33.  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
  34.  * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
  35.  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
  36.  */
  37.  
  38. #include <ansidecl.h>
  39. #include <math.h>
  40. #include <float.h>
  41. #include "ieee754.h"
  42.  
  43. /* Return the remainder of X/Y.  */
  44. __CONSTVALUE double
  45. DEFUN(__drem, (x, y),
  46.       double x AND double y)
  47. {
  48.   union ieee754_double ux, uy;
  49.  
  50.   ux.d = x;
  51.   uy.d = y;
  52. #define x ux.d
  53. #define    y uy.d
  54.  
  55.   uy.ieee.negative = 0;
  56.  
  57.   if (!__finite (x) || y == 0.0)
  58.     return NAN;
  59.   else if (__isnan (y))
  60.     return y;
  61.   else if (__isinf (y))
  62.     return x;
  63.   else if (uy.ieee.exponent <= 1)
  64.     {
  65.       /* Subnormal (or almost subnormal) Y value.  */
  66.       double b = __scalb (1.0, 54);
  67.       y *= b;
  68.       x = __drem (x, y);
  69.       x *= b;
  70.       return __drem (x, y) / b;
  71.     }
  72.   else if (y >= 1.7e308 / 2)
  73.     {
  74.       y /= 2;
  75.       x /= 2;
  76.       return __drem (x, y) * 2;
  77.     }
  78.   else
  79.     {
  80.       union ieee754_double a;
  81.       double b;
  82.       unsigned int negative = ux.ieee.negative;
  83.       a.d = y + y;
  84.       b = y / 2;
  85.       ux.ieee.negative = 0;
  86.       while (x > a.d)
  87.     {
  88.       unsigned short int k = ux.ieee.exponent - a.ieee.exponent;
  89.       union ieee754_double tmp;
  90.       tmp.d = a.d;
  91.       tmp.ieee.exponent += k;
  92.       if (x < tmp.d)
  93.         --tmp.ieee.exponent;
  94.       x -= tmp.d;
  95.     }
  96.       if (x > b)
  97.     {
  98.       x -= y;
  99.       if (x >= b)
  100.         x -= y;
  101.     }
  102.       ux.ieee.negative ^= negative;
  103.       return x;
  104.     }
  105. }
  106.